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Abstract. We present a new method to estimate time delays from light curves of lensed quasars. The method is 
based on minimization between the data and a numerical model light curve. A linear variation can be included 
in order to correct for slow long-term microlensing effects in one of the lensed images. An iterative version of 
the method can be applied in order to correct for higher order microlensing effects. The method is tested on 
simulated light curves. When higher order microlensing effects are present the time delay is best constrained 
with the iterative method. Analysis of a published data set for the lensed double Q 0957-1-561 yields results in 
agreement with other published estimates. 
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1. Introduction 

The time delay between light rays from gravitationally 
lensed quasar images is a measurable parameter directly 
related to the gravitat ional potential and to the Hubble 
constant Hq (Refsdal 1964 ). Accurate time delays ob- 
tained from multiply lensed QSOs can hence be used i) to 
determine Hq providing that the lens mass distribution is 
known, or ii) to constrain the mass distribution in a given 
lens, once Hq has been determined from other methods 
(or other lensed QSO systems). Much effort has therefore 
been devoted to observations of lensed QSOs during the 
last years, and in particular to the photometric monitoring 
of lensed QSO images. 

Measuring the time delay between the images of lensed 
quasars is, for several reasons, not a trivial task. First it 
requires regular monitoring of a target over a long period 
(substantially longer than the time delay). Second, the 
sampling is crucial and has to be determined on basis of 
the intrinsic variations in the quasar and the estimated 
time delay. Third, most objects are not observable during 
the whole year, i.e., they will be below the horizon for cer- 
tain periods. Fourth, there are nights with bad weather 
conditions when no data are obtained. Finally, variations 
in lensed quasars may not only be due to intrinsic fluctu- 
ations of the quasar, but also to microlensing by compact 
objects along the line of sight. Such a microlensing signal. 



depending on its timescale and amplitude, can be used 
to constrain the size of the continuum and the line emit- 
ting regions in the quasar, and the distribution of com- 



pact matter in the lens galaxy (Paczynski 1986, Kayser 
et al. 1986). However, as long as these external variations 
are not clearly distinguished from the intrinsic variations, 
microlensing remains a nuisance for time delay measure- 
ments. For the above reasons, advanced statistical meth- 
ods have to be used to measure time delays from quasar 
light curves. 

We have developed a method based on minimiza- 
tion of the data and a numerical model light curve. The 
aim was to develop a method based on simple principles, 
and to be able to properly detect and correct for possi- 
ble microlensing effects. It has already been successfully 
applied to several time delay measurements (Burud et 
al. 2000; Hjorth et al. [200 1| ). In this paper we will present 
the principles of the algorithm (Sect. ^ apply it to various 
simulated light curves (Sect. ^) and to a public dataset of 
QSO 09574-561 from Serra-Ricart et al. (p^99|) (Sect. |). 
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2. Method 

Let us assume that we have light curves for two images, 
A and B, of a lensed quasar. There are N data points 
in each of the time dependent light curves a{t) and b{t) 
with measurement errors (Ta and ab- The two curves are 
identical except for a shift in time, At, and in magnitude. 
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Am. It is therefore possible to model the two curves with 
one model curve g{t), and the parameters representing the 
time delay and the magnitude shift. In some cases the light 
from one or both of the images may be microlensed by 
individual stars in the lensing galaxy. Except for the cases 
of high amplification events, which are of short duration, 
microlensing effects are often slow variations that can be 
modeled to first order as a linear variation with slope a. 

As model curve we choose an arbitrary curve with a 
fixed number of equally spaced sampling points, M. This 
model may be minimized to the two observed light 
curves. The minimization is done only for the N observed 
data points. In this way only the model curve is interpo- 
lated and not the data. We thus minimize the following 
function: 



i=l 



^ 'ibiti-At)-{Am + ati))-g{tiy^ 



(1) 



The determination of the time delay from sampled 
light curves always relies on the implicit assumption that 
intrinsic variations are continuous in time and slow enough 
to be measured, given the adopted frequency of the obser- 
vations. Consequently with this assumption, if the typical 
sampling interval is n, we can thus smooth the model 
curve g{t) on the same time scale n (e.g., typically 7 days 
when observations are obtained once per week). This is 
done by introducing a smoothing term which minimizes 
the square of the difference between the original model and 
the same model convolved with a Gaussian r whose Full 
Width at Half Maximum (FWHM) is n. This smoothing 
is performed for all the points M on the model curve g(t) 
so that the dates without data points are smoothed out on 
the model curve and gaps in the observed light curves will 
not contribute significantly to the result. The smoothing 
term is multiplied by a Lagrange parameter A. This pa- 
rameter can be chosen so that the model curve matches 
the data correctly in a statistical sense (i.e., we get the 
correct x^). In practice however, excessive smoothing will 
prevent the minimization from converging. Therefore the 
smoothest possible solution which allows the minimization 
to converge is usually adopted. 

In a x^ minimization it is assumed that all the data 
points are independent. However, a light curve often con- 
sists in groups of nearby points and more isolated points. 
In regions where the time sampling is much shorter than 
the intrinsic quasar variations, the different measurements 
do not really bring new information on the shape of the 
light curve. Rather, they increase the precision with which 
the quasar magnitude is known at that moment. They 
might even be affected by common systematic errors, i.e. 
coming from a microlensing variation with a time scale 
comparable to the interval covered by these neighbouring 
points. 



In such cases of strongly uneven time coverage, we find 
that the best results are often obtained when a weight 
Wi is given to each data point in the following manner. 
For a point at a time U we calculate the relative distance 
ti — tj to all the other points in the curve. A Gaussian 
function with a FWHM= 2-\/2 ln(2)T2 chosen by the user 
is centered at U. The inverse of the sum of the ordinates of 
the Gaussian for each of the points tj will give the weight 
for the point U. The expression of Wi can be written: 



Wi 



1 



(2) 



The weight is normalized so that the maximum allowed 
weight for a data point is 1 , and this will occur only when 
one point is present within the time interval defined by 
T2. The ideal choice of T2 is the approximate time scale of 
variations in the quasar. For some systems this time scale 
will be long, e.g., several hundred days, for other systems 
it can be shorter e.g., 50 days. We can now write the final 
function to be minimized: 



N 



T = j2Wi 

i=l 

N 
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a{ti) - g{ti 



M 
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(3) 



where r * g signifies a convolution between r and g. 

2.1. Microlensing effects and iterative version of the 
algorithm 

Slow microlensing variations in one of the light curves are 
modeled as a linear term with the parameter a. If higher 

order variations are present we can use the method in an 
iterative way. This is done by splitting the light curves 
into several parts and analyzing them separately. The x^ 
method is applied to determine the amount of linear ex- 
ternal variations (a) for each separate part and for a range 
of input time delay values. First we run the programme 
with the first time delay dti in the chosen range and we 
estimate a in each of the separate parts. We then correct 
each part with the corresponding a and add the parts back 
together in order to obtain a set of microlensing-corrected 
light curves for the given input time delay dti. We can 
now run the x^ method to measure the time delay for 
these "corrected" curves. This procedure is repeated for 
all the input time delay values dt2, dt^, ... dtn in the range 
to be studied. We finally obtain a time delay measurement 
for each of the n input values, and all these mc^asurements 
will generally converge towards the real time delay value. 
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3. Time delay measurement on simulated light 
curves 

3.1. Data simulations 

We apply the method to simulated light curves in order 
to understand and estimate the errors on the measured 
values. Various real quasar light curves were used as mod- 
els for typical time variations and sampling of data points. 
The measurement errors in the curves represent simulated 
photon noise errors, so that the faint component will have 
larger errors than the bright one. We present the results 
from two sets of simulated light curves. The first set (la, 
b and c in Table |l|) consists of two curves, A and B, each 
of 60 points, dispersed over a time interval of 2 years (see 
Fig. |l| (top) and | (top)). The B curve is shifted by 145 
days in time and has a magnitude offset of 1.95 mag. In 
set lb we have added a slow linear time dependent vari- 
ation with a slope a (Eq. ^ in the B curve in order to 
simulate long-time-scale microlensing effects (Fig. |l| mid- 
dle). A higher order microlensing variation (at + j3t^) is 
added to the B curve in set Ic (Fig. |l| bottom). In the 
second set (2a, b, and c) we have simulated another kind 
of intrinsic variation and sampling. The curves span an 
interval of two years and they have 36 data points each. 
The time delay is 110 days and the magnitude offset is 
0.67 mag (see Fig. ^ top). In sets 2b and c, we added lin- 
ear (b) and higher order (c) variations in the A curve (see 
Fig. H middle and bottom). 

3.2. Time delay measurements 

The algorithm was applied to all the simulated sets of 
light curves. The fit was performed several times in order 
to check for the infiuence of various parameters such as: (i) 
the number of data points in the model curve (M in Eq. 
(it) the FWHM (t) of the Gaussian used to smooth the 
model curve {r{t) in Eq. |), (ni) the FWHM {2^2ln{2)T2) 
of the Gaussian that defines the time scale of the varia- 
tions in the curve (used to determine the weight Wi in 
Eq. and (iv) the Lagrange parameter (A in Eq. 
For the first set the best results were obtained when in- 
cluding a weight {Wi in Eq. ^ whereas for the second 
set the best results were obtained with no weight. This is 
because the data points in the second data set are more 
regularly distributed over time. The smoothing parame- 
ters were chosen so that the model curve is smooth over 
time, but still allowing variations large enough to fit the 
data. The results are shown in Table ^ 

We performed Monte Carlo simulations to estimate 
the errors in the results. Sets of 1000 light curves were 
simulated with error bars determined randomly from a 
Gaussian distribution with a tr equal to the measurement 
errors {a a and Ub). One set was made with the same sam- 
pling as the original curves (the simulated examples) and 
another set of 1000 curves was created with the same num- 
ber of sampling points as the original, but randomly dis- 
tributed. 
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Fig. 1. Top: The simulated light curves of two quasar com- 
ponents in set la (see text and Table [|). Middle: The time 
delay and magnitude shifted curves in set lb. The remain- 
ing difference between the two curves corresponds to the 
linear microlensing variation. Bottom: The light curves in 
set Ic, shifted in time delay and magnitude. The remain- 
ing difference between the two curves is the higher order 
microlensing variation. 
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Fig. 2. Top: The time delay and magnitude shifted Ught 
curves form set la. Bottom: The microlensing-corrected 
light curves derived by the iterative run on the curves in 
set Ic. 

When applying the iterative algorithm it is natural to 
think that the choice of the number of individual parts n 
is important for the results. Various tests on our simulated 
curves indicate that if n is too small, there is little differ- 
ence from the total curve and higher order microlensing 
effects can not be properly corrected for. If on the other 
hand n is too high, each part will contain very few data 
points and for some parts (depending on the sampling) 
the overlap between the two curves will be too small to 
determine the external variations. The choice of n will thus 
depend on the number of data points and on how they are 
distributed (e.g., presence of gaps in the curves). The opti- 
mal values for the curves in the simulated curves, i.e., the 
values resulting in the best convergence of At, were found 
to be 3 and 2 in sets 1 and 2 respectively. The results 
shown in Table |l| correspond to the mean and standard 
deviation of the values obtained. 
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Fig. 3. Top: The simulated light curves of two quasar com- 
ponents in set 2a (see text and Table [|) . Middle: The time 
delay and magnitude shifted curves in set 2b. The remain- 
ing difference between the two curves corresponds to the 
linear microlensing variation. Bottom: The light curves in 
set 2c, shifted in time delay and magnitude. The remain- 
ing difference between the two curves is the higher order 
microlensing variation. 
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Fig. 4. The R-band light curves of QSO 0957+561. The 
curve is shifted by 423 days and 0.0625 mag. The solid 
line shows the model curve derived by the ^ algorithm. 

3.3. Results 

All the time delay measurements obtained in sets la, b 
and 2a, b are very well determined within the estimated 
errors (see Table Q). As expected, the results depend on 
the sampling of the light curves and the error estimates 
from the Monte Carlo simulations with random sampling 
are larger than the ones obtained with a fixed sampling. 
In sets lb and 2b where a slow linear variation is added in 
one of the curves, the time delay is still well determined 
whereas in sets Ic and 2c where higher order variations are 
added there are small systematic errors in the time delay 
estimated with the direct method. In these cases, the it- 
erative method gives more accurate results with smaller 
uncertainties. Fig. ^ displays how the microlensing effects 
in data set Ic are corrected for with the iterative method. 
In the sets with no microlensing (la and 2a) we note 
that the iterative method yields larger errors than the di- 
rect method. When no external variations (simulated mi- 
crolensing effects) are present the best constrained time 
delay value is obtained by setting a = 0, hence removing 
one of the free parameters. The iterative method should 
therefore preferably only be applied where external varia- 
tions are clearly present. 

The magnitude shifts between the two curves are very 
well determined in sets 2a, b and c whereas in example 
la, b and c there are small systematic errors of 0.01 — 0.05 
mag. 

4. Application to light curves of QSO 0957+561 

Our method was also applied to a public data set o f 
QSO 0957-^561 published by Serra-Ricart et al. ( |l999|) 
(see Fig. ^). The R-band lightcurves contain 214 data 



points and span the time interval from February 1996 
to July 1998. We removed eight points that deviate by 
more than 0.05 mag from neighbouring points in both 
components. We applied our algorithm to the data 
and, as was done for the simulated data, we estimated 
the error bars using Monte Carlo simulations. A time de- 
lay estimate At = 423 ± 9 days and a magnitude shift 
Am = 0.063 ± 0.007 were obtained. Considering the two 
gaps in the data set, JD 2450242 to 2450347 and JD 
2450637 to 2450729, Serra-Ricart et al. have divided the 
data set into two parts, DS I and DS II. Since the sec- 
ond part (DS II) show the largest activity they have de- 
termined a time delay value on DS II only, correspond- 
ing to the autumn 1997/spring 1998 season. They mea- 
sure a time delay value of 425±4 days and a magnitude 
shift of 0.06 mag on these curves. Measurements with two 
other methods, the dispersion spectra and discrete cross- 
correlation techniques give values of 426 ± 12 and 428 ± 8 
days respectively on the same data set (Serra-Ricart et al. 
1999). In order to compare the results we therefore also 
apply our method to DS II. On the curves in DS II we 
measure At — 424 ± 22. Both time delay values, the one 
measured on the entire data set and the one measured on 
DS II, are compatible with the results from Serra-Ricart 
et al. However our errors are larger than the ones found 
by Serra-Ricart et al. In particular we note that our er- 
rors are considerably larger when using only DS II than 
when using the entire data set. Serra-Ricart et al. do not 
give any time delay measurement from the total curve but 
their method seems to be well suited for time delay mea- 
surements on a continuous curve (DS II), even when the 
time span of observation is only of the order of the time 
delay value. With our method the time delay is better con- 
strained with longer time series, even if gaps are present 
in the curves. 

An additional reason for the differences in the results 
might be that whereas we have discarded 8 points form 
the total light curve, Serra-Ricart et al. have discarded 23 
points. Given our limited knowledge of the data points we 
found no reason to exclude more points. 

5. Conclusion 

We have developed a method to determine time delays 
between light rays from lensed quasar images from their 
respective light curves. As other based methods (e.g.. 
Press et al. 1992 and Barkanna 1997) we assume that the 
statistics of the measurement errors are Gaussian, and we 
require smoothness for nearby points in the light curves. 
In our method we have included an optional weight to each 
of the data points, depending on the relative distance to 
neighbouring points. In this way one point isolated in time 
will receive more weight than each individual point in a 
cluster of nearby points. 

Since microlensing effects are often encountered in real 
light curves of lensed quasars we have adapted our method 
in order to correct for such variations. A linear term in one 
of the curves is included to correct for slow long-term vari- 
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Table 1. Results from the analysis of the simulated light curves. The sets (a, b and c) correspond to simulations 
with respectively no, linear and higher order external variations in one of the components. These variations simulate 
microlensing effects (see the text for more details). 



light curves 


la 


lb 


Ic 


2a 


2b 


2c 


true 


145 


145 


145 


110 


110 


110 


derived At: 














direct method 


143±5 


151±5 


132±5 


109±3 


119±5 


119±4 


M.C. fixed sampling 


146±5 


148±5 


135±5 


109±3 


114±5 


115±4 


M.C. random sampling 


155±23 


149±16 


147±18 


118±35 


117±32 


118±24 


Iterative method 


145±7 


145±2 


143±9 


108±6 


111±3 


114±8 


true Am 


1.950 


1.950 


1.950 


0.670 


0.670 


0.670 


derived Am; 














direct method 


1.923±0.002 


1.968±0.004 


1.868±0.004 


0.667±0.001 


0.667±0.003 


0.647±0.004 


M.C. fixed sampling 


1.927±0.002 


1.964±0.004 


1.885±0.004 


0.667±0.001 


0.665±0.003 


0.678±0.004 


M.C. random sampling 


1.921±0.005 


1.963±0.036 


1.905±0.08 


0.672±0.017 


0.65±0.04 


0.65±0.02 


Iterative method 


1.90-2.0 


1.91-2.01 


1.95-2.31 


0.64-0.67 


0.65-0.69 


0.65-0.69 



ations. Realistic simulations show that the time delay and 
the magnitude shift between two light curves are deter- 
mined to within a few percent for cases with no or slow 
microlensing effects. For faster microlensing variations we 
find that running the algorithm in an iterative way yields 
better time delays. This confirms what was proposed in 
the time delay measurement of B1600-I-434 (Burud et al. 



200C ) for which four methods were applied: the minimum 
dispersion method (Pelt |l996| ), the SOLA method (Pijpers 



1994, 1997 ), the method described in the present paper 
and its iterative version. The iterative method yielded the 
best constrained time delay in this data set due to external 
variations in one of the components. 

We also applied our method to a public data set of 
QSO 0957-1-561 and have shown that our results are in 
agreement with the published time delay. 
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